function Plot_Bond_Risk_Premia(parms,gesol,Phi_Nxz,prob_Nxz)

    % select relevant range
    P1 = 0.001; P2 = 0.999;
    idx_N1 = Inf; idx_N2 = -Inf;
    idx_x1 = Inf; idx_x2 = -Inf;
    
    for iz = 1:parms.Nz
        cdf_N = cumsum(sum(prob_Nxz(:,:,iz),2)); cdf_N = cdf_N/cdf_N(end);
        idx_N1 = min(idx_N1, find(cdf_N<=P1,1,'last'));
        idx_N2 = max(idx_N2, find(cdf_N>=P2,1,'first'));
        cdf_x = cumsum(sum(prob_Nxz(:,:,iz),1)); cdf_x = cdf_x/cdf_x(end);
        idx_x1 = min(idx_x1, find(cdf_x<=P1,1,'last'));
        idx_x2 = max(idx_x2, find(cdf_x>=P2,1,'first'));
    end
    

    rp_2y = calc_rp(gesol,Phi_Nxz,24);
    rp_5y = calc_rp(gesol,Phi_Nxz,60);
    
    figure; m = 2; n = 2; idxplot = 0;
    
    phi_of_x = @(x)1./(1+exp(-x));
    [NN,PP] = meshgrid(parms.grid_N(:),phi_of_x(parms.grid_x(:)));
    
    for iz = 1:parms.Nz
        idxplot = idxplot + 1; subplot(m,n,idxplot);
        mesh(NN,PP,rp_2y(:,:,iz)')
        xlabel('N')
        ylabel('\phi')
        zlabel(['rp, 2y, z=',num2str(iz)])
        set(gca,'xlim',parms.grid_N([idx_N1 idx_N2]))
        set(gca,'ylim',phi_of_x(parms.grid_x([idx_x1 idx_x2])))
        
        idxplot = idxplot + 1; subplot(m,n,idxplot);
        mesh(NN,PP,rp_5y(:,:,iz)')
        xlabel('N')
        ylabel('\phi')
        zlabel(['rp, 5y, z=',num2str(iz)])
        set(gca,'xlim',parms.grid_N([idx_N1 idx_N2]))
        set(gca,'ylim',phi_of_x(parms.grid_x([idx_x1 idx_x2])))
    end

end

function rp = calc_rp(gesol,Phi_Nxz,mat)
    
    idx_buy = find(gesol.real_bonds.maturities==mat);
    idx_sell = find(gesol.real_bonds.maturities==mat-12);
    idx_1y = find(gesol.real_bonds.maturities==12);
    
    p_buy = log(gesol.real_bonds.prices(:,:,:,idx_buy));
    size_p = size(p_buy);
    p_sell = log(gesol.real_bonds.prices(:,:,:,idx_sell));
    p_1y = log(gesol.real_bonds.prices(:,:,:,idx_1y));
    
    rp = reshape((Phi_Nxz^12)*p_sell(:),size_p) - p_buy + p_1y;

end